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Modeling individual mobility’s impact on COVID-19 transmission 


Abstract 


This research explores the influence of individual mobility on COVID-19 
transmission, utilizing a temporal mathematical model to clarify disease 
spread and vaccination dynamics across diverse regions. Employing a com- 
putationally efficient two-patch configuration that emphasizes regional in- 
teractions, our study aims to guide optimal disease control strategies. The 
introduced SEIR-V model with a two-patch setup estimates the vaccination 
reproduction number, R,, while equilibrium points and system stability 
are identified. Visualizations from numerical simulations and sensitivity 
analyses illustrate key parameters affecting the vaccination reproduction 
number and COVID-19 control measures. Our findings underscore system 
responsiveness, emphasizing the intricate relationship between Ry, migra- 


tion rates, and disease prevalence. 


AMS subject classifications (2020): Primary 39A12; Secondary 92C60, 92D30. 


Keywords: COVID-19; Metapopulation; Global stability; Local stability; 


Vaccination reproduction number. 


1 Introduction 


In December 2019, a coronavirus emerged in Wuhan, China, transmitted from 
animals to humans [22]. The outbreak escalated rapidly, with 44 reported 
cases in China by January 3, 2020 [29]. The virus swiftly spread nationwide, 
prompting the declaration of a state of epidemic in March 2020, exacerbated 
by global human migration and travel between cities [18, 20]. The facilitation 
of worldwide travel and commerce contributed to the rapid global dissemi- 
nation, resulting in over 585 million reported cases and 6.4 million reported 
deaths as of August 2022 [28]. Consequently, a multitude of researchers have 
been captivated by the dynamics of COVID-19 dissemination [16, 6, 7, 9, 12]. 

The investigation focused on the movement patterns of individuals across 
a network during the global spread of COVID-19 in 203 countries [13]. 
Researchers observed that human mobility, particularly through cities and 
tourism, significantly contributed to the virus’s dissemination, supported by 


Arino and others who affirmed that epidemic spread increases during trans- 
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portation [2]. Numerous epidemiologists, including Baister et al. [4], Anmad 
et al. [1], McCarthy et al. [21], and Iyaniwura et al. [15] employed metapop- 
ulation models to explore COVID-19 transmission dynamics. These models 
considered diverse scenarios, such as studying individual mobility and social 
interactions in major cities, investigating vaccination strategies in a prison 
setting, and using hybrid gravity-metapopulation modeling to analyze the 
spread between different regions. Other studies delved into specific regions 
like Ireland and formulated generalized n-patch SEIR epidemiological models 


to devise effective control strategies against the disease [14, 19]. 


This paper introduces metapopulation modeling to investigate the impact 
of individual mobility on the spread of COVID-19, building upon our previ- 
ous model [7]. Our focus is on developing a temporal mathematical model to 
describe COVID-19 transmission with vaccination, aiming to optimize disease 
control strategies within and between different regions. Using a two-patch 
configuration (n = 2), we address the complexities of calculations, ensuring 
the boundedness and positivity of solutions. We calculate the vaccination re- 
production number (R,) for the two-patch metapopulation and analyze the 
local and global stability of the disease-free equilibrium (DFE), establishing 
that if R, > 1, then the disease becomes endemic. Sensitivity analysis and 
numerical simulations further explore the impact of key parameters on Rv, 
confirming the existence and stability of the endemic equilibrium. Addition- 
ally, we examine the influence of migration rates on COVID-19 transmission 


between the two patches. 


The paper’s structure unfolds as follows: Section 2 revisits the previously 
introduced and analyzed model. In Section 3, we introduce the metapopula- 
tion model. Section 4 proposes the SEIR-V model with a two-patch configu- 
ration for COVID-19, estimating the vaccination reproduction number (R,). 
We identify equilibrium points and conduct stability analysis for the dynamic 
system of COVID-19 propagation. Section 5 presents numerical simulations 
and sensitivity analyses to illustrate theoretical results. Finally, Section 6 


encapsulates our conclusions. 
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2 Baseline model 


We first recall the model introduced and analyzed in [7], which describes the 
transmission dynamics of COVID-19. In this model, the population is sub- 
divided into seven distinct groups: susceptible individuals (S), exposed but 
not yet infectious individuals (£), vaccinated individuals (V), asymptomatic 
infectious individuals (I,), symptomatic and hospitalized infectious individ- 
uals (I,), individuals who have succumbed to the disease (D), and those who 
have recovered (R). 

It is assumed that there is no vertical transmission of the disease so that 
all births occur into the susceptible human class, at the rate 4: > 0. Moreover, 
we assume that the total human population N remains constant; thus, birth 
and death rates are equal. 


Hence, the normalized reduced system is given by 


S!(t) = w— (aan + a2 (1 — )) S(t) Ua (t) + Is (t)) — WS (Et) 
+AgV (t) + yR(t) — ALS (t) , 
Et (t) = (ain + a (1 — 9) S(t) Ua (t) + Ts (t)) — (81 + B2 ++ rz) E(t), 
V' (t) = ALS (t) + ABE (t) + Aaa (t) — An +A5 +H) V (it), 
Ty, (t) = BE (t) — (11 + Aa t+ 4) a (4), 
Ty (t) = BoE (t) — (v2 + 6+) Is (t), 
R' (t) = ya (t) + Yals (t) + AsV (t) — yR(t) — pR(E), 
(1) 
where 
a is the infection rate of confined susceptible and az is the infection rate 
of unconfined susceptible such as a, < a2. 
7 is the confinement rate within the population. 
(3;);-1,9 are infected-infectious rates. 
y is the reinfection rate after recovery from a first infection. 
(Yi);=1,2 are the recovery rates. 
(Ai) ;=1,3,4 are the vaccination rates. 
Ag is the rate of vaccine ineffectiveness. 


As is the vaccine effectiveness rate. 
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6 is the COVID-19 mortality rate. 
i is the natality rate. 
The study of the model is detailed in [7]. For this model (1), the vacci- 


nation reproduction number is given by 


(ay + a2 (1 —1)) 
Tes Kaa) Gant Dae) ©) 
(Bi (V2 +6 + Bh) + Bo (M1 FAG + B)) Ao t+ As +H) (7+ E) 


(81 + Bo + wt 3) ((A2 +A + 1) (¥ +H) + AL (7 +A5 +H) 


3 Metapopulation model 


In works [4, 1, 21, 15, 14, 2, 19, 3], a comprehensive system of differential 
equations is developed to characterize human mobility. Within this model, 
subpopulations are identified based on their origin and their current location. 
Building upon the model detailed in the preceding section, we expand it 
by introducing the concept of a neighborhood wherein interactions occur 


between infected and non-infected individuals. 


Consider a network comprising p nodes. Within this network, individuals 
are characterized by two key attributes: the node of origin, signifying their 
residence, and the node they currently occupy at time ¢. We assume that 
the total human population N; in each node remains constant and strictly 
positive, and the transition rates from one pathological state to another also 
remain the same for all patches. Furthermore, we suppose that birth occurs 
in the resident node while deaths take place in any node where the human is 
present. 

The population in the patch 2, is divided into compartments of suscepti- 
ble (.S;), exposed (;), vaccinated(V;), asymptomatic infectious (I,,;), symp- 
tomatic and hospitalized infectious (I;,;), death (D;), and the recovered as 
(Ri). 


The total number of individuals in the patch 7 size at time t is defined by 


N; = S;(t) + Ei(t) + Vi(t) + Ia i(t) + Io,i(t) + Ri(t) + Di(t), 


and the total population is given by 
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We adopt a methodology akin to that employed by Arino and van Den 
Driessche [3]. Migration is observed between any pair of patches, transpiring 


at the specified rates: 


e pe migration rate of susceptible individuals from the patch 7 to the 


patch j, 


° pk, migration rate of infected individuals (noninfectious) from the patch 


i to the patch J, 


° Yh; migration rate of vaccination individuals from the patch 7 to the 


patch j, 


° gis migration rate of asymptomatic infectious individuals from the 


patch i to the patch J, 


° yg) j migration rate of symptomatic infectious individuals from the patch 


i to the patch J, 


° ph migration rate of recovered and immunized infectious individuals 


from the patch 7 to the patch J, 


where 


8 E Vv Te Ts R 
Pi = Pu = Pu = Pi Pii pi, = 0. 


Hence, the metapopulation system of differential equations that will model 


the dynamics of coronavirus spread is the system given by 
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dS; Pp 
BEM ~ (aim to2(1—m)) Tait La) Sit > PhS; 
J=1(9 4%) 
Pp 
_ ( > é) Si — (wt A1) Si + AV; + 7Ri, 
J=1(9#%) 
di => (aN; + ag (1 = ™)) (Lai = I, 4) S; “fF » iL; 
j=1G#4i) 
p 
-( eB) B—(.4 tne 
j=1(GFi) 
dV; P 
ran MS; + ASE; +Mlait SO vf; 
j=1(G#i) 
Pp 
< ( PS é) Vi-— (Ao +A5 +H) Vi, 
j=1(GFi) 
dla, P I 7 I 
aie AB + YS gittag - Dy Pie | Lae — (1 + Aa t+ HB) Las, 
jJ=1 (Gi) j=l (G#i%) 
dls j 
oar BoE; — (72 +04 p) Is, 
dR; P 
a Vitae + Yats,4 + ASV + 2 pk R; 
J=1 (Gi) 
Pp 
ze ( > é) Re +4) Ry 
j=1G 4%) 
dD; 
: = ols im) 
dt , 


with initial conditions 
S; (0) > 0, E;, (0) > 0, yy E;,(0) > 0, 
V; (0) > 0, Iq; (0) = 0, Pen Iq; (0) > 0, 
Is 3 (0) = 0, baie Is4 (0) > 0, R; (0) 2 0, 


for i =1,p. 


4 Two-patch SEIR-V model for COVID-19 


To enhance the realism of our study, we confine our investigation to a two- 


patch model, undertaking a rigorous mathematical analysis of its dynamics. 
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The diagram of the transmission of COVID-19 in two-patch is shown in 
Figure 1. 


Figure 1: Diagram of the proposed two-Patch SEIR-V Model. 


The SEIR-V two-Patch model is formulated as follows: 
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dS 
- = pny — C1 Jaa + Is1) Si + y3, Se = (v2o ela co) Si + AM + Ri, 
dS S Ss 
“dt = png — C2 (La,2 + Is.2) So+ pi251 = (v3 as co) So + AgV2 + yRo, 
dE 
aE =¢) (La + 1,1) Si + ph Ee — (yp) +3) Ei, 
dE. 
= = 9 (Lao + Is,2) S2 + PEL — (yh, +3) Eo, 
dV, Vv Vv 
he = A181 + A3E, aie Nala or p31 V2 _ (Yl2 + ca) v1, 
dVo V Vv 
ae Ai S29 + Ag Ee + Agla,2 + PioVi — (en + c) Vo, 
dIq 
= 6 F, + glttna — (els + cs) leds 
dla 
aE = Pio + oily — (ek + cs) La,2, 
dl, 
ae = Bok) — cols, 
dl, 
i = BoE — cols.2, 
dR 
a" = Yaa + yale +AsVi + ph Ro — (ib + e7) Ri, 
dR 
a = W112 + Yets.2 + A5sV2 + yi Ri a (os ag C7) Ro, 
(3) 
with initial conditions 
Ej, (0) > 0, Fy (0) + Fy (0) > 0, 
Lai (0) S 0, Tq 1 (0) + Ta2 (0) > 0, 
Is (0) Z 0, I, 1 (0) =F Iso (0) 2 0, 


for i= 1,2. Here 
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Ni No 

n ng = ’ 

1 ee N ca = (Ag +A5 4+ 1), 

co = ’ 

0 i cH = trAat+p, 

a = (aim +a2(1—m)), C6 = 7 o+ pu oy 
6= 72+ 

c2 = (ain2 + a2 (1—72)), c7 = (y+) 
7= + : 


c3 = By + Bo +pt+Asz, 


4.1 Boundedness and positivity of solutions 


We consider the following zone of biological interest: 
= { ($1, 52, Bi, Basi, Vas Tats Ta,2s Tots Toa» Ray Re) € RY? : (6) 
O< So 5+ Bt OVit+ So lait Do doe t+ D> Ri < rental, 


Theorem 1. Consider system (3) with initial conditions (4). 


The set 2 is an attracting and positively invariant with respect to model 
(3), with $;,V;,R; (i= 1,2); remaining positive. The total population in 
each patch is assumed to be constant, and all their pathological states are 


also bounded. 
Proof. Under initial conditions (4): 


e In the event that E, becomes zero at t; before Ey becomes zero, then 
at t1, we have 


dE 
a = ¢1 (Jaa + Jo,1) 51 +o, Es > 0. 
This demonstrates that EF, is a nondecreasing function of t at t;. Hence, 


EE, stays nonnegative. Analogously, the same holds for Fo. 


e In the event that [,,1 becomes zero at some time tz before Ig,2 becomes 
zero, then at tz, we have 


dla1 
dt 


=F, + yatta, > 0. 


This demonstrates that J, is a nondecreasing function of t at to. 


Hence, J,,1 stays nonnegative. Analogously, the same holds for Iq,2. 
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e In the event that J, ; becomes zero at some time t3 before J,.2 becomes 


zero, then at ts, we have 


dl. 1 
—_ = BoE, > 0. 
a Bo Ey 


This demonstrates that J, is a nondecreasing function of ¢ at ts. 


Hence, J,,; stays nonnegative. Analogously, the same holds for J, 9. 


e Suppose now that at some time t4, S; (t4) = 0 before Sp goes to zero. 


Then at t = t4, from system (3), we have 


$2 + A2V, + 7R1 > 0, 


as 


dt 
Thus, there is no time t, such that $ (t4) = 0. Therefore, 5 stays pos- 


N. 
which implies that > 0 when ra is strictly positive. 


itive for t > 0 when the initial condition S; (0) > 0. Using comparable 


reasoning, we deduce the positivity of So. 


e Likewise at some time ts, Vi (ts) = 0 before V2 goes to zero. Then at t 


= ts, from system (3), we have 


dV; 
a = 1S) + A3E1 + Nala + Yd V2 > 0, 


d 
which implies that = > 0 because S; (t) > 0 for t > 0. 


Thus, there is no time ts such that V; (ts) = 0. Therefore, Vi stays 


positive. Using comparable reasoning, we deduce the positivity of V2. 


e Likewise at some time te; Ri (te) before Rz goes to zero. Then at t¢ 


= te, from system (3), we have 


dR 
_ = 11Ta1 + y2le1 +AsVi + oh Re > 0, 


which implies that “a > 0 because Vj (t) > 0 for t > 0. 


Thus, there is no time tg such that R; (tg) = 0. Therefore, Ri stays 


positive. Using comparable reasoning, we deduce the positivity of Ro. 
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Since the positive set Q is invariant under system (3) and the total population 


in each patch (N;),;_, 5 is assumed to be constant, and all their pathological 


states are also bounded. 


4.2 Disease-Free Equilibrium 


The DFE (disease-free equilibrium) of system (3) is expressed as 
i 
Ej = (St Sp00Vs VF 0000Ri Rs) , 


where 


gt = pkoks (ny (d, + dg) + dine) 
- dg (dz + dy + da) 
gt = pkoks (dani + (dz + d3) n2) 
? d3 (d3 + d, + da) 
At (ph + c4) Lkgk3 (nq (dy + d3) + dinz) 
kad3 (dg + dy + dz) 
4 Parka ikaks (dyn, + (dz + dg) nz) 
ds (d3 + d, + do) : 


= A vlopkeks (ny (dy, + ds) + ding) 
7 kodz (d3 + dy + dz) 
Mt (C26 + ca) bkgk3 (dgn1 + (dz + d3) nz) 
k2d3 (d3 + dy + dz) ; 
A5Arkapkoks (ny (dy, + d3) + din2) 
kok3d3 (d3 + d; + dz) 


Vy — 


VS 


Ri = 


Lkok3A5\1 (psi ks + psiiks) (don1 + (do + ds) nz) 


kok3d3 (yi, + e7) (ds + di + de) 
Rt = Lk2k3X15 (yioks + pth ka) (ny (dy - d3) + dyn) 


kok3d3 (pet, + C7) (ds + dy Tr dz) 
AsArkspkaks (dani + (dz + d3) n2) 
kok3d3 (d3 + d; + dz) 


dy = yA\1Aske + AA2kseyy + kgkophh, 
dy = y\Askz + \Anks yo + kaksypo, 
d3 = k3kou+ Ap (A5 + C7) (ke +k7+ c4C7) ‘ 
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and 
2 Vv 
kg =cg+ YYoCa Tr Yo1C4, 


kg = 7 + pier + pacer, 

ka = (oh +. c4) (95, + 7) + ohh, (9) 
ks = (plo + ca) (91 + e7) + oben, 

ke = yn (934 +7) + 084 (vio +4), 


kr = pp (pis + cr) + 9B (on + c4) - 


The existence of the DFE €5 remains unaffected by any prerequisites regard- 


ing positive operational data. 


4.3 Vaccination reproduction number R, of two-patch 


metapopulation 


In this subsection, we present the method used for our SEIR-V model to esti- 
mate the vaccination reproduction number R,. This method is proposed by 
Diekmann, Heesterbeek, and Metz [26] and elaborated by Van Den Driessche 
and Watmough [27, 25], which gives a way of determining R, for an ordi- 
nary differential equation compartmental model by using the next-generation 


matrix. 


Let X = (Fi, Eo, La, La,2, 15,1, La Then the system can be written as 


dX 
OX = F(x) - V(X), 
where 
e1 Jaa + Is) $1 
c2 (La,2 + Is,2) So 
0 
F(X) = ‘ 
0 
0 
and 
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E E 
—~p Bg + (vi ae c3) Ey 


— pF, + (pH + ¢3) Eo 
—B, Ey — tlie p13 C5 ) La 


V(X)= 
—Bi Ee — pi3laa + (p33 +65 | Ine 
— BoE + cel s1 
—B2E2 + cels.2 


The basic reproduction number, R,, is calculated by next generation tech- 


nique. The F' and V matrices at the DFE € is given as follows: 


00 cS} 0 cSt 0 
00 O C255 0 C255 


_ @F(X)| _[00 0 0 0 0 
— OX lee 100 0 0 0 0 |’ 
00 0 0 0 0 
000 0 0 0 
and 

ph+c —pe, 0 0 00 
-yf, pita 0 C OO 
vy V(X)| _| -f& 0 vate —~H 0 0 
OX ex 0 —B1 £13 pai tes 0 0 
—Be 0 0 0 0 
0 =; 0 0 06 


la I 
* * e1(cs+¢54 ) * C1P2i Ox C1 Ox 
miSf m3S7 ST i) ee ae 


ky 
maS3miS3 Beiss SCF) ss 9 ass 
Ae Oo. 2 0 0 0 0 | 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
where 
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ere (‘< + 3) (= a =A) % Pu 


C6 ky c3ky1 


bt vg, (Bo. csBi\ , Bivis 
ko C6 ky c3ky ; 


E I 
_ ya (P2 C51 Arpat 
calla ( ko (2 ky ) " c3ky4 , 


(cz + Ya) (2 =P) | at) 


ko C6 ky e3ky 


and 


a E E 2 
ko = C3Pj9 + c3pa1 + €3, 


ky = cspi3 + espa} + cB. 
The reproduction number, R,, is the spectral radius of the next generation 
matrix, which is given by 


_ 1 Lkoks 
~ 2d3 (d3 +d, + dz) 


+y((a + n d3) my, — (do + n2d3) ma)? + Amoms (dy + n,d3) (do + nat) ; 


Ry 


(Ca + nid3) my, + (do + ngd3) ma 


where (c;) and (d;) are defined in (5) and (8), respectively. 


4.4 Model analysis 


Theorem 2. If R, < 1, then the DFE €j is locally asymptotically stable 


and unstable otherwise. 


Proof. System (3) has a unique DFE point €j in the set ©, given in (7). 


The eigenvalues of the Jacobian matrix at the DFE point €} are given by 


Pil = P2 = —C6,; 
03 = —(%ia + oh +02), 
pa = —C3, 

= Ta Ta 
ps =— (cs + 9194+ oii) , 
Po = —C5, 
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where (/;);—7 12 are roots of 
PYZ)\= 2° £0,294 C17* 40,7" 46,277 £0, 2 + Cy. 


By using Lienard—Chipart criteria, the DFE point € is locally asymptotically 
stable if the coefficients C5, C4, C3, Co, Ci, and Co are positive and their 
Hurwitz’s determinants are positive. It shows that the positivity of Hurwitz 
determinants is very complicated to achieve. However, the computational 
cost can be reduced by using the simplified version mentioned in [5, 17]. 
That is, the coefficients that must be computed are C5, Cy, C3, Co, C1, and 
Co only. Note that Cs, C4, and C3 are needed because the simplified version 
of Lienard—Chipart criteria still require |H2| = C5C4 — C3 > 0 [11]. 


Indeed, we employed formal calculation software (Maxima Version 5.42.1) 
to establish that if R, <1, then all (C;);95 and |Hp| are positive. 


Theorem 3. If R, < 1, then the the DFE € is globally asymptotically 
stable in 2. 


Proof. Since for 1 = 1,2, we have 


51 < ST, S2 < 53, 


dE 
= < = (ph + c3) Ey + ph, E> + cSt Lat + a St 1s, 
dE: 
_s < pi F — (ph +3) Bo + c2S$Ia2 + c2S$Is,2, 
dIq, 

< Ak, — (vis + cs) Ta, + 9341a,2; 

(10) 

dl g,2 Ta a 

dt < Bi 24+ vi5la1 — | Goi +5) La2; 
dIs4 

; < E 3 I, , 

dt Be 1— C6Lls1 
dIs.9 
— < BoE — cels.o. 

ae Be 2— Céts,2 


Defining an auxiliary linear system using the right-hand side of system (10), 


we have 
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om = — (ef +03) E+ ph Ea t+ a Silas + Sits, 
ure = phi — (phi + cs) By + coS$3la2 + e253 Is,2, 
ee = BEY = (vis vr cs) Ta + erila,2, 

a = B, By + giglaa — (ek + cs) Ia, 

7 = BoE, — celei, 

- = BoE — cele, 


or, in other words, 


dels. ass es a, ys. 
= (& Eg La La,2 15,1 2) 


dt 
—(yiet+e3) ph cS% 0 ast 0 Ey 
re = (phi + c3) 0 C255 0 255 Ey 
_ By 0 _ (vis + cs) Owi3 0 0 ia 
- 0 by els -(eh+e) 0 0 |] hes 
Bo 0 0 0 —cg 0 Tac 
0 Bo 0 0 0 —c¢% To 
So, 7 7 
dE, dE, 
dE» dE» 
ON ey | Se (11) 
dt dao dIy 2 
dIs.1 dl. 
dI.,2 di, 2. 


We have R, < 1 <= > a(F —V) <0, where o (F — V) is the spectral abscissa 
of the matrix F' — V (see the proof of Theorem 2 [26]). So, when R, < 1, 
the eigenvalues of (F — V) are with negative real parts. Thus all nonnegative 
solutions of (11) are such that 


lim (A, Ep, jeer Tas T..1, i.2) = (0, 0, 0, 0, 0, 0) : 


t+ 
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By a standard comparison principle (see [24, Theorem B.1]) and the non- 
negativity of (£1, Es, Ia1,1e,2,1s1,Is,2), we conclude that when R, < 1, all 
nonnegative solutions of (3) satisfy 


lim (£1, Eo, Lai, Ia,2, 15,1, [s,2) = (0, 0,0,0,0,0). (12) 


t—>+00 


Since (12) is satisfied, (3) is an asymptotically autonomous system [see [8, 


Theorem 2.5]] with limit affine system 


ds: 
a = pt v3\S2 — (Po +c0) Si + %2Ni + 7Ri, 
dS: 
oa = + Yi2S1 — (3, + co) S2 + A2V2 + 7Ro, 
dV, 
“y= Sa + ehiVe — (va + ca) Mis 
(13) 
se er y 
a 192 + Yio — (oh as ca) Va, 
dR 
=a =AsVi + ph Ro — (yp +7) Ri, 
dR 
a = AsVo + vi, Ri — (vs) + €7) Ro. 


It is easy to verify from (13) that 


_iim ($1, 52,Vi, V2, Ri, Ro) = (ST, 53,Vj Va , RT, Rp) : 


The globally asymptotic stability of the DFE €% is then proved if R, < 1. 


The presence of endemic equilibrium, characterized by positive infective 
numbers, has barely been discussed here. Although their existence has not 
been formally proven, numerical simulations suggest a globally asymptoti- 
cally stable equilibrium when R, > 1. Investigating the system’s persistence 
properties under the condition R, > 1 would be a valuable avenue for further 


exploration. 


Lemma 1. [30] Let ¢; : X — X be a semiflow and let Xo C X be an open 
set. Define 0X9 = X \ Xo and Mg = {x € OX: da € OXo,t > 0}. Assume 


the following conditions: 


1. d:Xo C Xo and ¢; has a global attractor A. 
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2. There exists a finite sequence M = {M,..., M;,} of disjoint, compact, 


and isolated invariant sets in 0X9 such that 


(a) Q (Mo) := Usemyw (x) C UE Mi. 

(b) No subset of M forms a cycle in 0X0. 

(c) M; is isolated in X. 

(d) Ws (M;)N Xo = 0, where W* (M;) = {x € Xo: w(x) C Mj} for 
each 1 <i<k. 

Then ¢; is uniformly persistent with respect to (Xo, 0X0), that is, there exists 


8 > 0 such that 
lim inf d(¢i2%,0X0) > 0 for x € Xo. 
t—>-+00 


Therefore, we have the following result. 
Theorem 4. If R, > 1, then system (3) is uniformly persistent, namely, 
there exists a constant # > 0 such that 


lim inf S;(¢)>6, lim infF£;(t)>6, lim infV;(t) > 89, 
t—+ +00 t—>+00 


t—+0o 


lim infl,;,(¢)>0, lim infl,;(t)>0, lim infR;(t) >, 
t—+ +00 ‘ t—+ +00 t—++0o0 


for any initial conditions satisfying (4). 
Proof. Choose 
x= 


Xo = {(S1, S2, 1, Be, Vi, Va, La,1, 10,2, 1,1, Is,2, Ri, Ro) € X 
Ey + Ey >0,1a1 + Ia,2 > 0,1s,1 + Is,2 > 0}, 


and 


OXo =X~ Xo 
= { ($1552, Bi, Ba, Vi, Vas Tats Ta,2s To,15To,2» Ri, Re) Ee Xx 


f= Beh =e Sa = 0}. 


Let ¢; be the semi flow induced by the solutions of system (3). 
We have proved in Theorem / that ¢:X9 C Xo and ¢; is ultimately 


bounded in Xo; so there always exists a global attractor for ¢,. It is obvious 


Tran. J. Numer. Anal. Optim., Vol. 14, No. 2, 2024, pp 580-612 


599 Modeling individual mobility’s impact on COVID-19 transmission 


that €9 is the unique boundary equilibrium on 0X0, which implies that €> 
is globally stable on OXo. Moreover, ($1, 52,,Vi, V2, Ri, Ro) converges to 
(97 Sasa Vis Vay 7, R3) on OX. 

Let M, = {Ej} and M = {Mj}, then Uzeu,w (x) = M; and no subset 
of M forms a cycle in 0Xo. 

If Ry > 1, then €§ is unstable in Xo. Therefore, conditions (2 — c) and 
(2—d) of Lemma 1 are satisfied. Then ¢; is uniformly persistent with respect 
to (Xo,0Xo0), that is, there exists 6 > 0 such that 


lim inf d(¢.%,0Xo0) 20 for x € Xo. 
t—+0o 


We have thus shown that if R, > 1, then the disease is endemic. 


5 Numerical simulation 


In this section, we extend our exploration into the numerical simulation 
of COVID-19 transmission within a two-patch framework, leveraging well- 
established methodologies commonly employed in epidemiological studies. 
Building upon the model detailed and examined in Section 4, our objective 
is to unravel the intricate dynamics of COVID-19 transmission influenced by 


migration, vaccination, and medical unavailability. 


In a prior study [7], we investigated disease transmission dynamics in the 
absence of migration in a single patch. In the current work, which builds 
upon our earlier research, we now consider the pivotal role of migration— 
representing the movement of individuals between patches. Migration rates, 
denoted by yi2 and 21, exert a substantial influence on the overall trans- 


mission dynamics. 


To ground our simulations in real-world scenarios, we integrate data from 
Table 1 into our analysis. Additionally, we set migration rates from the 
patch 1 to the patch 2 (y?, 9%, vo, pis, yt) at 0.001 and migration 
rates from the patch 2 to the patch 1 (v1, 944, 0%, ps, pe) at 0.005. 
These rates underscore the interconnectedness of patches and contribute to 


the comprehensive understanding of disease spread across different regions. 
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Table 1: The operational parameters of the simulation consist of two categories: those 


characterizing the population and the disease remain the same across all patches, 


whereas parameters reflecting the behavior of individuals vary from one patch to another 


Parameters Description ecrae 
N Population size 100000 
Ni Population size in patch 1 (Most populated patch). 75000 
N2 Population size in patch 2 (Least populated patch). 25000 
ie Recruitment rate. 3.38¢-05 
a1 Infection rate for unconfined persons. 0.01-0.3 
a2 Infection rate for confined persons. 0.01-0.3 
"1 The confinement rate in patch 1. 0.2-0.5 
12 The confinement rate in patch 2. 0.2-0.5 
Bi Rate of an exposed becoming asymptomatic infected. 0.092 
po Rate of an exposed becoming symptomatic infected. 0.05 
Ay Vaccination rate of susceptible. 0.01-0.03 
Ao The ineffectiveness rate of the vaccine 0.001765 
A3 Vaccination rate of exposed persons. 0.00167 
A4 Vaccination rate of asymptomatic infected. 0.0167 
As The effectiveness rate of the vaccine. 0.008 

y Rate of a recovered becoming susceptible. 0.0099 
y1 The recovery rate of asymptomatic people 0.01 

yo the recovery rate of symptomatic people 0.0005 
O COVID-19 death rate. 0.0055 


Our simulation parameters encompass scenarios without migration re- 


strictions, incorporating measures such as containment and vaccination strate- 


gies. The outcomes of these simulations gauged through R, and other perti- 


nent metrics, offer valuable insights into the potential efficacy of interventions 


aimed at curtailing the virus’s propagation. 


Subsequently, we delve into an analysis of the sensitivity of key parame- 


ters to the vaccine reproduction number R,. Following this, we meticulously 


scrutinize the theoretical results previously obtained, aiming to provide em- 


pirical substantiation and validate the accuracy of our theoretical framework. 
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This multifaceted approach enhances our understanding of the practical im- 
plications of proposed interventions and reinforces the robustness of our the- 


oretical foundation. 


5.1 Sensitivity analysis of the main parameters on R, 


To effectively mitigate the impact of COVID-19 and reduce associated fatali- 
ties, understanding the relative importance of various parameters influencing 
its progression is crucial. The vaccination reproduction number, R, [10], 
stands out as a key determinant in the transmission dynamics of COVID-19. 

Our investigation employs sensitivity analysis on R, within the framework 
of our model (3). This analysis aims to discern the individual importance of 
each parameter, shedding light on influential factors in disease transmission. 
The derived sensitivity indices offer insights into how a state variable changes 
relative to variations in specific parameters. 

Sensitivity analysis, a widely used technique, involves computing the first- 
order partial derivatives of model output concerning input parameters. This 
process can be conceptualized as calculating gradients in a multidimensional 
reference parameter space [23], providing valuable information on the model’s 


responsiveness to slight changes in input parameters. 


The normalized forward sensitivity index of R,,, denoted as ine quantifies 
the sensitivity with respect to a given parameter p and is defined by the 


expression: 


pe 
P Op R, 


This index facilitates a quantitative assessment of how alterations in individ- 


ual parameters influence the value of R,, a pivotal metric for understanding 
disease transmission. By leveraging explicit formulas for each parameter af- 
fecting R,, as detailed in reference [31], we gain precise insights into the im- 
pact of specific parameters on R,. Such sensitivity analysis proves invaluable 
in prioritizing interventions and control measures to effectively curtail the 
spread of the disease, enabling a more targeted and comprehensive approach 
to managing COVID-19. 
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Table 2: Local sensitivity indices for R, concerning parameters in the proposed model are 
computed at the baseline parameter values provided in Table 1. Parameters exhibiting 
the highest sensitivity are highlighted in bold 


Parameters a, Qe 1, Ny Bi Bz A 

Sensitivity Index 0.11 0.89 -0.11 -0.00006 -0.34 0.37 -0.18 

Parameters Az A3 Ag As Y V1 Y2 

Sensitivity Index 0.032 -0.016 -0.018 0.066 0.079 -0.11 -0.059 
Iq 

Parameters P21 P21 P21 P24 P21 


Sensitivity Index 0.15 0.0003 0.008 0.002 0.009 


Parameters Pi2 Piz Pie 915 Oe 
Sensitivity Index -0.15 -0.006 -0.008 -0.009 -0.009 


Table 2 highlights the parameters most responsive to changes in the vac- 
cination reproduction number (R,), with the transmission probability of 
unconfined individuals (a2) and the transition rate from exposed to symp- 
tomatic infected (82) exhibiting the highest sensitivity. The vaccination rate 
(Ai) and the transition rate from exposed to asymptomatic infected (1) 
closely follow in significance. A 10% increase (or decrease) in a2 or (2 leads 
to corresponding 8.9% and 3.7% changes in R,, respectively, as denoted by 
TR = +0.89 and T» = +0.37. Conversely, a 10% increase in d, or fy 
corresponds to a 1.8% and 3.4% decrease in R,, respectively, indicated by 
TX” = —0.18 and Tj" = —0.34. 


The local analysis emphasizes the importance of implementing effective 
containment measures for managing the transmission probability (a2) and 
the transition rate from exposed to symptomatic infected. This involves 
stringent containment measures, the use of physical barriers, and adherence 
to social distancing guidelines. Additionally, vaccination efforts, represented 
by the rate (A1), play a crucial role in reducing the spread of the disease. By 
strategically controlling these parameters, a significant reduction in COVID- 
19-related deaths can be achieved. 
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5.2 The stability of the DFE 


In this subsection, we will confirm the theoretical results obtained earlier 


through numerical simulations for the following parameter values (per day): 
a, = 0.010 ; ag = 0.0655 ; 7, = 0.20 ; ng = 0.25 ; Ay = 0.01817. 
The existence and stability conditions in Theorem 2 are satisfied 
Ry = 0.8375 < 1. 


Hence, system (3) has a DFE, 
* * * * * * * * * * * ok * vi 
& = [St S3, Ey, Ey, VW; Wy, Ta Ty 21 Lea» ee Ri, Rj] 
= (0.1915, 0.0385, 0, 0, 0.3326, 0.0713, 0, 0, 0, 0, 0.2861, 0.0573)" . 


Therefore, the DFE € 9 is stable. This also implies that the disease is eradi- 
cated, as depicted in Figures 2 and 3. 

Furthermore, if we increase the values of the most sensitive parame- 
ters such as ay = 0.0653, ag = 0.1295, 4, = 0.20, m2 = 0.30 and 
A, = 0.00117, through direct calculation, the vaccination reproduction num- 
ber R, = 6.4251 > 1. Consequently, according to Theorem 2, the DFE is 


unstable, and the model (3) has an endemic equilibrium. 


Patch 1 Patch 2 


———Susceptible 
Vaccinated 
Recovered 


Uninfected states 
=) 
& 


0 i 0 i i 
0 1000 2000 3000 4000 5000 0 1000 2000 3000 4000 5000 
Time (in Day) Time (in Day) 


Figure 2: The proposed model has only one DFE &€, which is asymptotically stable 
when Ry, < 1. 


Tran. J. Numer. Anal. Optim., Vol. 14, No. 2, 2024, pp 580-612 


Bouziane, Boubekeur, Keddar and Belhamiti 


Infected states 


604 


x10° Patch 1 x10" Patch 2 
4 
— Exposed 
—— Asymptomatic Infectious 
i —'—'*Symptomatic and Hospitalized Infectious P 
2 
i 
4 
{ 
h : 
1 i MWe seddbe bed ONDA eS ed eee Mba tended baad Aahe Heine aN dele we we 
iJ 
\ 
X: 
: , : : 
1000 2000 3000 4000 5000 0 1000 2000 3000 4000 


0 
0 


Time (in Day) 


5000 
Time (in Day) 


Figure 3: The proposed model has only one DFE €%, which is asymptotically stable 


when Ry < 1. 


5.3 The existence and stability of the endemic 


equilibrium of the model 


Now, we deal with the existence of endemic equilibrium of model (3). Con- 


sider the following parameter values (per day): 


Hence, we obtain 


System (3) has an endemic equilibrium, 


a 4 # # # # # 
Ei* = [Si*, So*, Ei*, Eo", Vi* , Vo™, 


Figures 4 and 5 illustrate the stability of €[*. 


Ry, = 6.4251 > 1. 


[** 


a,l> 


ek 
Ty, 


a1 = 0.0653, ag = 0.1295, m1 = 0.20, no = 0.30, A, = 0.00117. 


T 
4 # # # 
Isis hey Ry ’ Ry | 


= [1.0663 — O01, 2.3693e — 02, 6.0569e — 04, 6.4518e — 06, 


1.6194e — 02, 3.0613e — 03, 2.0238e — 03, 8.2520e — 05, 
T 
5.0185e — 03, 5.3556e — 05, 1.5170e — 02, 2.7128e — 03] : 
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Patch 1 Patch 2 
0.8 7 7 7 : 10.25 : 
: : : -——Susceptible | i : : 
OT here : See ee : Se eee —Vaccinated f 
" : ——Recovered | 9.9 H Pccmansinaee hclate round Wenaentnutios wintnend wnmenuciss | 
06 t FETC TETOT RTT PLOT ENT TTT Ter Tea tLI Tee PIERO 1 
\ 1 
3 1 \ : P 
5 05 : tet etrti crt rete ere ety ete 4 0.15 La. semis : ome : sult pera ms m| 
\ 4 
ie} I! 1 : 
2 WER ccercte Sy apennsanrernastrapemsanciarglecperesamnseueayets 1 
2 it \ 
S 01 oak Keenan. eerrenrenst rte 
ti 03 oe es Tre Per rte 
pb ! \ 
1 
0.05 bes 


5 ‘ 
0 1000 2000 3000 4000 5000 0 1000 2000 


3000 4000 5000 
Time (in Day) 


Time (in Day) 


Figure 4: Illustration of the stability of €y*. 


Patch 1 a 


x10 Patch 2 
0.18 H i T I eee TESST ESE ET] 
BY —— Exposed j : H : 
a6 iN “ —— Asymptomatic Infectious | 
o14be 4s —'='* Symptomatic and Hospitalized Infecti 


Infected states 


: ; 6 males 
0 1000 2000 3000 4000 5000 0 1000 2000 3000 
Time (in Day) 


4000 5000 
Time (in Day) 


Figure 5: Illustrates the stability of E]*. 


By adjusting the key operational parameters to yield varying values of R, 
greater than 1, we have 


Ry, = (3.7816, 4.6628, 5.5440, 6.4251, 7.3063, 8.1875)’ . 


We observe in Figure 6, the endemic equilibrium point along with standard 
deviations for each disease state, 
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SD= [3.6400 — 02, 7.1634e — 03, 3.8114e — 05, 4.1860e — 07, 
4.1299e — 03, 8.3947e — 04, 1.2731e — 04, 5.2295e — 06, 
T 
3.1569e — 04, 3.4790e — 06, 3.1924e — 03, 6.5984e — 04] . 


These standard deviation values for each disease state indicate that the data 
points in a set are closely clustered around the mean, suggesting low vari- 


ability. 


Additionally, numerical simulations in Figure 6 show that the solutions 
with different values of R, converge to the same endemic equilibrium, where 


the initial values are 


S; (0) = 7.1970e — 01, $2 (0) = 2.3990e — 01, E; (0) = 8.2500e — 04, 
Ep (0) = 2.7500e — 04, V; (0) = 2.7975e — 02, Vo (0) = 9.3250e — 03, 
Tq,1 (0) = 3.7500e — 04, Iq.2 (0) = 1.2500e — 04, I,.1 (0) = 3.7500e — 04, 
I,.2 (0) = 1.2500e — 04, Ry (0) = 4.5000e — 04, Ro (0) = 1.5000e — 04. 
O14 T T T T T 


+ +k +k +k *k Kk #K +k +e +k #k #K 
Sp Sy EY BY Vy Va Tay a2 Ig g2 Ry 


Figure 6: Standard deviation values for each disease state for E]*. 
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Certainly, there is a possibility that the endemic equilibrium remains 
stable given the specified parameters. This leads us to an intriguing question: 
if the vaccination reproduction number R, exceeds 1, then the model (3) with 


migration allows for a locally asymptotically stable endemic equilibrium. 


5.4 The effect of migration rate on the transmission of 
COVID-19 in the two patches 


For migration rate values ranging from 0 to 0.01, we observe the results in 


Figure 7. 


Vaccination reproduction number R,, 
0.01 


0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 
The migration between patch 1 and patch 2 


oa oO Oo oo oa 
[=] i=) f=] i=) 2 
oS i=) =) oa a 
Wa n ~~ co Bd 


Oo 
So 
Ss 
& 


Oo oO 
i=] So 
oa oa 
ie) tes} 


The migration between patch 2 and patch 1 


vi 


Figure 7: The migration rate effect on Ry. 


The numerical simulations presented in Figure 7 underscore the endemic 
nature of the disease. As the migration rate diminishes from the more popu- 
lated patch 1 to the less populated patch 2, and concurrently rises from the 
patch 2 to the patch 1, the vaccine reproduction number Rv demonstrates 
a proportional increase. This observation validates the outcomes of the sen- 


sitivity analysis performed on Rv. Theoretical insights, coupled with these 
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numerical simulations, emphasize the intricate relationship between COVID- 
19 prevalence across diverse patches and the corresponding migration rates. 
These findings contribute to a comprehensive understanding of the dynamics, 
suggesting that the prevalence is intricately tied to the patterns of migration 


between distinct patches. 


6 Conclusion 


In conclusion, the findings from sensitivity analysis shed light on the nuanced 
interplay between key parameters and the vaccination reproduction number 
(R,) in the context of COVID-19 dynamics. Notably, the transmission rate 
of unconfined individuals (a2) and the transition rate from exposed to symp- 
tomatic infected individuals (G2) emerge as pivotal factors influencing the 
increase in Ry. Robust containment measures targeting these aspects, such 
as stringent protocols, physical barriers, and adherence to social distancing 
guidelines, are imperative for effective disease control. Furthermore, the vac- 
cination rate (A) and the transition rate from exposed to asymptomatic 
infected individuals (61) closely follow in significance. A 10% increase or 
decrease in a2 or G2 corresponds to substantial changes in R,, emphasizing 
the sensitivity of the system. Conversely, a 10% increase in A, or 8; leads to 
a notable decrease in R,. These results provide actionable insights into the 
impact of parameter variations on the vaccination reproduction number. 
For adequate operational data, the stability of DFE €3 suggests the poten- 
tial for disease eradication. However, elevating the values of sensitive param- 
eters leads to an increase in R, to 6.4251, surpassing the critical threshold 
of 1. This instability of €) indicates the presence of an endemic equilib- 
rium, as per Theorem 2. The theoretical framework aligns with numerical 
simulations, as depicted in Figures 2 and 3, underscoring the importance of 
both analytical and computational approaches in understanding the system 
behavior. The low standard deviation values across different disease states 
suggest that the data points are closely clustered around the mean, indicat- 
ing limited variability. Additionally, numerical simulations in Figures 4 and 
5 further affirm the endemic nature, demonstrating that varying values of 


R, > 1 converge to the same endemic equilibrium. 
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We conclude the simulation by highlighting the complex relationship be- 


tween COVID-19 prevalence across distinct patches and the corresponding 


migration rates, offering valuable insights into the dynamic nature of the 


disease and guiding strategies for effective disease management. 
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